      program main
c
c     This is an example program to test the symmetric skyline solver.
c
c     Input:
c            n		dimension of the problem
c            nb		average band width desired when generating data
c            nv         the std. dev. wanted for the band width
c            igen       = 0      the matrix is generated randomly
c                       > or = 1 the position of the diagonal elements
c                                must be given

      implicit double precision (a-h,o-z)

      integer          id(10001)
      double precision a(1000000),ao(1000000),b(10000)

100   print *, ' '
      print *, 'Enter n, average bandwidth, bandwidth variance:'     
      read( 5, *, end=999 ) n, nb, nv
      print *, 'Storage Scheme: Enter  1  for NORMAL Jennings '
      print *, '                       2  for REVERSE Jennings:'
      read( 5, *, end=999 ) ityp
      print *, 'Matrix to be generated randomly: '
      print *, '       Enter  0  for whole matrix generated randomly'
      print *, ' '
      print *, '              1  for diagonal elements given,' 
      print *, '                     matrix generated randomly'
      print *, ' '
      print *, '              2  for diagonal elements given,'
      print *, '                     matrix element given :'
      read( 5, *, end=999 ) igen

c     Use as many CPU as available
      ncpu = mp_numthreads() 

c     Use Jennings Method because reverse Jennings is NOT
c     available yet.
      print *, ' '
      if( ityp .EQ. 1 ) then
	  print *,'main:  jenning\'s method '
      else
	  print *,'main:  reverse jenning\'s method '
      end if

      print *,'       num cpu in system = num cpu used: ', ncpu

      call driver( n, nb, nv, igen, ityp, a, ao, b, id )

      go to 100

999   stop
      end

      subroutine driver( n, nb, nv, igen, ityp, a, ao, b, id )

      implicit double precision (a-h,o-z)

      integer           n, nb, nv, igen, ityp
      integer           id(*)
      double precision  a(*),ao(*),b(*)

      integer           am
      integer           iops
      real              etime,t1,et,wt,tarray(2),et1,second
      real*8            rand
      double precision  tol
      double precision  zero, one, two, half, onem
      double precision  tmp

      parameter         ( zero = 0.0d0 )
      parameter         ( one  = 1.0d0 )
      parameter         ( two  = 2.0d0 )
      parameter         ( half = 5.0d-1 )
      parameter         ( oneh = 1.0d2 )
      parameter         ( onem = 1.0d6 )

      integer           idmach
      double precision  dmach
 
      intsize = idmach( 'I' )
      iflsize = idmach( 'F' )
      tol     = dmach( 'E' )

      if(igen .eq. 0) then
        iseed=10001
        call srand(iseed)
        do i=1,n
           a(i) = rand()
        enddo

        id(1)=1
c	if( ityp .EQ. 1 ) then
           do i=2,n
              col= abs(two*(a(i)-half)*float(nv)+float(nb)+half)
              ilen=min0(int(col),i)
              ilen=max0(ilen,1)
              id(i)=id(i-1)+ilen
           enddo
c	else
c	   id(2) = 2
c           do i=2,n
c              col= abs(two*(a(i)-half)*float(nv)+float(nb)+half)
c              ilen=min0(int(col),i)
c              ilen=max0(ilen,1)
c              id(i+1)=id(i)+ilen+1
c           enddo
c        end if

      elseif(igen .ge. 1) then
           read(5,*) (id(i),i=1,n)
      endif

      am=id(1)
      hm=one/float(id(1))
      mx=id(1)

      do i=2,n
        j=id(i)-id(i-1)
        mx=max(j,mx)
        am=am+j
        hm=hm+one/float(j)
      enddo

      print *, ' '
      print*,'       ndof: ',n
      ip=100.*float(mx)/float(n)+.5
      print*,'       max band: ',mx,' or ',ip,' % '

      am=am/float(n)
      iam=am
      ip=100.*am/float(n)+.5
      print*,'       arithmetic mean band: ',iam,' or ',ip,' % '

      hm=hm/float(n)
      hm=one/hm
      ihm=hm
      ip=100.*hm/float(n)+.5
      print*,'         harmonic mean band: ',ihm,' or ',ip,' % '

      lena=id(n)
      i=float(8*lena)/onem+.5
      print*,'       num of terms: ',lena,' megabytes: ',i

      if(igen .le. 1) then
        do i=1,lena
           a(i) = rand()
        enddo
        do i=1,n
          a(id(i))=(a(id(i))+half)*oneh
        enddo
      elseif(igen .eq. 2) then
        read(5,*) (a(i),i=1,lena)
      endif

      do i=1,lena
        ao(i)=a(i)
      enddo
      do i=1,n
        b(i)=one
      enddo

      if( ityp .EQ. 2 ) then
	  id( n+1 ) = id( n ) + 1
	  do j = n, 2, -1
	     ibgn = id( j-1 ) + 1
	     iend = id( j )
	     id( j ) = ibgn
	     imax = ( iend - ibgn + 1 )/2
	     do i = 1, imax
		tmp = a( iend )
		a( iend ) = a( ibgn )
		a( ibgn ) = tmp
		ibgn = ibgn + 1
		iend = iend - 1
	     end do
	  end do
      end if

         et1=second()
         t1=etime(tarray)
         call dskydc(a,n,id,info,ityp)
         if(info .gt. 0)  then
	    print *, '**** INFO = ', info
	    return
	 end if
	 if(info .LT. 0) return

         et=etime(tarray) - t1
         wt=second()-et1



      call djesol(a,n,id,b)
      call error(ao,n,id,b)

      big=zero
      do i=1,n
        err=abs(b(i)-one)
        if(err .gt. big) then
          big=err
        endif
      enddo
      print*,'       error= ',big
      
      if(nv .eq. 0) then
        xb=nb
        xn=n
        ops=xb*(xb-1.)*(3.*xn-2.*xb+1.)/3.
        ops=ops/onem+.5
      else
        iops=0
        call counter(id,n,iops)
        ops=float(iops)/onem+.5
      endif
      if(et .gt. zero) flops=ops/et
      if(wt .gt. zero) then
         wflops=ops/wt
         ratio=wt/et
      endif

      write(6,48) ops
48    format(/,' million fp ops = ',f8.2)
      write(6,49) et,wt
49    format(' cpu time  = ',f6.2,' wall time  = ',f6.2)
      write(6,50) flops,wflops,ratio
50    format(' cpu mflops= ',f6.2,' wall mflops= ',f6.2,' ratio= ',f6.2)

      return
      end

      SUBROUTINE COUNTER(LC,NEQ,ICOUNT)
C---------------------------------------------- 
C      ---ACTIVE COLUMN EQUATION SOLVER  AX=B-- 
C          KK=1 LDL FACTORAZATION ONLY
C          KK=2 FORWARD REDUCTIION ONLY 
C          KK=3 BACKSUBSTITION ONLY 
C---------------------------------------------- 
      DIMENSION LC(NEQ)
      DOUBLE PRECISION S,T,FF
      ICOUNT = 0
 50   IF(NEQ.EQ.1) RETURN 
C------LDL FACTORAZATION--------------------- 
      DO 500 J=2,NEQ
      JH=LC(J)-LC(J-1)
      IF(JH.EQ.1) GO TO 500 
      K=J-JH+1
C-----FORM U(I,J) - TOP OF COLUMN DOWN TO DIAGONAL
      I=K 
      LCIM1=0
      IF(I-1 .GT. 0) LCIM1=LC(I-1)
 100  NT=MIN0(JH-J+I,LC(I)-LCIM1)-1 
      IF(NT) 300,300,150
 150  NS=LC(I)-NT 
      NE=LC(I)-1
      IJ=LC(J)-J+I
      IC=IJ-LC(I) 
C**      S=0.0D0 
      IF(I.EQ.J) GO TO 400
      ICOUNT = ICOUNT + 2 * NT 
C**      CALL DOTP(A(NS),A(NS+IC),S,NT)
C**      A(IJ)=A(IJ)-S 
 300  I=I+1 
      GO TO 100 
C-----FORM L(I,J) AND U(I,I) ------------------ 
 400  DO 450 N=NS,NE
      ND=LC(K)
      K=K+1 
C**      T=A(N)
C**      A(N)=A(N)/A(ND) 
      ICOUNT = ICOUNT + 3
C** 450  S=S+A(N)*T
 450  CONTINUE
C**      A(IJ)=A(IJ)-S 
      ICOUNT = ICOUNT + 1
 500  CONTINUE
      RETURN
      END 
 
      SUBROUTINE ERROR(A,N,ID,X)
      POINTER (PLHT,LHT), (PZ,Z)
      INTEGER  N,ID(1),LHT(1),M
      DOUBLE PRECISION A(1),X(1),Z(1),S,TOL
      INTEGER          IDMACH
      DOUBLE PRECISION DMACH

      INTSIZE = IDMACH( 'I' )
      IFLSIZE = IDMACH( 'F' )
      TOL     = DMACH( 'E' )

      PLHT=MALLOC(N*INTSIZE)
      PZ=MALLOC(N*IFLSIZE)
      LHT(1)=0
      DO 5 I=2,N
        LHT(I)=ID(I)-ID(I-1)-1
    5 CONTINUE
      DO 10 K=1,N
        Z(K)=X(K)*A(ID(K))
   10 CONTINUE
      DO 40 K=2,N
        M=ID(K-1)+1
        II=0
        S=Z(K)
        DO 30 I=K-LHT(K),K-1
          Z(I)=Z(I)+X(K)*A(M+II)
          S=S+X(I)*A(M+II)
          II=II+1
   30   CONTINUE
        Z(K)=S
   40 CONTINUE
      DO 50 I=1,N
        X(I) = Z(I)
   50 CONTINUE
      CALL FREE(PLHT)
      CALL FREE(PZ)
      RETURN
      END

